 #include <math.h>
#include <iostream>
#include "dexelfield.hpp"

const double NaN=0;

long total;

dexelfield::dexelfield()
{
  dx=0;dy=0;
  xmin=0;
  xmax=0;
  ymin=0;  
  ymax=0;
  zmin=0;
  zmax=0;
  
  nx=0;
}

dexelfield::dexelfield(double _xmin, double _ymin,
		       long _nx, long _ny,
		       double _dx, double _dy)
{  
  dx=_dx;
  dy=_dy;
  xmin=_xmin;
  ymin=_ymin;

  zmax=-1000;
  zmin=1000;

  nx=_nx;
  
  xmax=xmin+nx*dx;
  ymax=ymin+_ny*dy;

}


void dexelfield::blank(double xmin, double xmax,
	double ymin, double ymax,
	double _zmin, double _zmax,
	double dx, double dy ){
 
 long mx=(xmax-xmin)/dx+1;
 long my=(ymax-ymin)/dy+1;

 *this=dexelfield(xmin,ymin,mx,my,dx,dy);

 zmin=_zmin;
 zmax=_zmax;
 
 for(long i=0; i<=mx*my; i++)
   f.push_back(ourdexel(_zmin,_zmax));

  total=f.size();
 
}

void dexelfield::adddexel(long ix,long iy,ourdexel d)
{
  f[hash(ix,iy)]=d;

  if(zmin>d.min())
    zmin=d.min();
  if(zmax<d.max())
    zmax=d.max();
}

long dexelfield::hash(long ix,long iy) const
{
  return iy*nx+ix;
}

long dexelfield::referx(double x) const
{
  return (long)round((x-xmin)/dx);
}

long dexelfield::refery(double y) const
{
  return (long)round((y-ymin)/dy);
}

double dexelfield::unrefx(long ix)
{
  return double(ix)*dx+xmin;
}

double dexelfield::unrefy(long iy)
{
  return double(iy)*dy+ymin;
}

long dexelfield::hash(double x, double y) const
{
  return hash(referx(x),refery(y));
}

body_state dexelfield::state(double x, double y, double z) 
{
  if( x<xmin || x>xmax ||
      y<ymin || y>ymax )
    return AIR;
  size_t has=hash(x,y);

  if(has>=total)
    return AIR;  //////////////!!!!!!!!!!!!!!!!!!!!!!FIXFIXFIX
  
  return f[has].state(z);
}

body_state dexelfield::state(long ix,long iy, double z) 
{
  if( ix<0 || ix>(xmax-ymin)/dx ||
      iy<0 || iy>(ymax-ymin)/dy )
    return AIR;
  size_t has=hash(ix,iy);

  if(has>=total)
    return AIR;  //////////////!!!!!!!!!!!!!!!!!!!!!!FIXFIXFIX
  
  return f[has].state(z);
}


double dexelfield::intersect_x(double x0, double y0,double z0,
			       double x1) 
{


  long iy0=refery(y0);
  long ix0=referx(x0);
  long ix1=referx(x1);
  
  if(ix0==ix1) return (x0+x1)/2;
  
  if(state(x0,y0,z0)==state(x1,y0,z0))
    return NaN;
  
  
  long ix=0;
  
  while(labs(ix1-ix0)>1)
    {
      body_state state0=state(ix0,iy0,z0);
     
      ix=ix0+(ix1-ix0)/2;
            
      body_state cstate=state(ix,iy0,z0);
      
      if(cstate!=state0)
	ix1=ix;
      else
	ix0=ix;
    };
   
  return xmin+dx*ix;
}

   
double dexelfield::intersect_y(double x0, double y0,double z0,
			       double y1)
{
 
  long ix0=referx(x0);
  long iy0=refery(y0);
  long iy1=refery(y1);

  if(iy0==iy1) return (y0+y1)/2;
  if(state(x0,y0,z0)==state(x0,y1,z0))
    return NaN;
  
  long iy=0;
  
  while(labs(iy1-iy0)>1)
    {
      body_state state0=state(ix0,iy0,z0);
	
      

      iy=iy0+(iy1-iy0)/2;
            
      body_state cstate=state(ix0,iy,z0);

      
      if(cstate!=state0)
	iy1=iy;
      else
	iy0=iy;
    };
   
  return ymin+dy*iy;
}

/*double dexelfield::updatezmin()
{
  zmin=0;
  
  typedef map<long,dexel>::iterator P;
  for(P p=f.begin(); p!=f.end(); p++)
    if(p->second.min()<zmin)
      // zmin=p->second().min();
  return zmin;
}

double dexelfield::updatezmax()
{
  zmax =0;
  
  typedef map<long,dexel>::iterator P;
  for(P p=f.begin(); p!=f.end(); p++)
    if(p->second.max()>zmax)
      // zmax=p->second().max();
  return zmax;
  
}
*/
  
double dexelfield::intersect_z(double x0, double y0,double z0,
			       double z1)
{

  if( x0<xmin || x0>xmax ||
      y0<ymin || y0>ymax )
    return NAN;
  size_t has=hash(x0,y0);
  
  if(has>=total) return NAN;
  
  ourdexel &d=f[has];
  
  // if(d.state(z0)==d.state(z1))
  //   return NaN;

  return d.border(z0,z1);
  
 
}


void dexelfield::apply_cutter(cutter *c)
{
  if(c->xmin()>=xmax || c->xmax()<=xmin ||
     c->ymin()>=ymax || c->ymax()<=ymin)
    return;
  
  
  double x_min=(xmin>c->xmin())?xmin:c->xmin();
  double x_max=(xmax<c->xmax())?xmax:c->xmax();

  double y_min=(ymin>c->ymin())?ymin:c->ymin();
  double y_max=(ymax<c->ymax())?ymax:c->ymax();

  double pc=0;
  int percents=0;
  double pp=
    (referx(x_max)-referx(x_min))*
    (refery(y_max)-refery(y_min));
  
  
  for(long ix=referx(x_min); ix<=referx(x_max); ix++)
    for(long iy=refery(y_min); iy<=refery(y_max); iy++)
      {
	interval cut=c->intersect_z(unrefx(ix),unrefy(iy));

	size_t has=hash(ix,iy);
	
	if(has<total)
	  f[has].substract(cut);

	 pc++;
	  
      }
   
}

	
	
	

  
      
